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ABSTRACT 

Millisecond pulsars (MSPs) are known as highly stable celestial clocks. Nevertheless, 
recent studies have revealed the unstable nature of their integrated pulse profiles, 
which may limit the achievable pulsar timing precision. In this paper, we present a 
case study on the pulse profile variability of PSR J1022+1001. We have detected ap¬ 
proximately 14,000 sub-pulses (components of single pulses) in 35-hr long observations, 
mostly located at the trailing component of the integrated profile. Their flux densities 
and fractional polarisation suggest that they represent the bright end of the energy 
distribution in ordinary emission mode and are not giant pulses. The occurrence of 
sub-pulses from the leading and trailing components of the integrated profile is shown 
to be correlated. For sub-pulses from the latter, a preferred pulse width of approxi¬ 
mately 0.25 ms has been found. Using simultaneous observations from the Effelsberg 
100-rn telescope and the Westerbork Synthesis Radio Telescope, we have found that 
the integrated profile varies on a timescale of a few tens of minutes. We show that 
improper polarisation calibration and diffractive scintillation cannot be the sole reason 
for the observed instability. In addition, we demonstrate that timing residuals gener¬ 
ated from averages of the detected sub-pulses are dominated by phase jitter, and place 
an upper limit of ~ 700 ns for jitter noise based on continuous 1-min integrations. 
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1 INTRODUCTION 


'ed pc__ 

lestial clocks (e.g. lYerbiest et al.ll2009l ) and thus excellent 
tools in the ongoing end eavours of high precision gravity 


tools m the ongoing end eavours ot nign precision gravity 
tests ( Hramer_et al. 200 (f), neutron sta r equations of state 
llDemorest et al.l 2010| ; [OzeLgLal] [2010|) , and gravitational 
wave detection llvan Haasteren et al. ffeoill: lYardlev et al.l 


l201ll : ll)emorest et al.ll2013l : IShannon et alj|2013h . These ex¬ 

periments utilise the high precision timing of millisecond 
pulsars (MSPs), which are known for their rotational stabil¬ 
ity and the stability of their integrated pulse profiles (which 
do not change over many years). However, recent studies 
have shown that the profiles of MSPs can exhibit low-level 
instabilities that may limit the achievable timing precision. 
In general, these instabilities can be categorized into two 
types: temporally uncorrelated variabilities caused by the 


instability of single pulses (e.g . Ijenet et~aklIl998l : iLiu et all 
120121 : IShannon fe Cordesl 12012 ), and tempor ally correlated 
varia t ions induced by other phen omena (e.g. iKramer et al.l 
1 19991 : [Edwards fc Stappersll2003al ). The contribution of jit¬ 
ter to the pulsar timing noise could be mitigated by simply 
increasing the exposure time or by correc ting the biased inte¬ 
grated profiles with suitabl e algorithms (ICordes fe Shannonl 
|2010| ; lOslowski et al.l 120131 ) . Comprehensive studies of pro¬ 
file instabilities would give new insights into the fundamental 
timing limits and thereby potentially improving the timing 
precision, which is essential for the aforementioned astro- 
physical experiments. 

PSR 11022+1001 is a MSP with a rotation period of ap¬ 
proximately 16 ms. It exhibits a regular rotational behaviour 
and is included in the current pulsar timing array observa- 
tions for the pur pose of gravitational wave detection (e.g. 
lManchestein2013l ). The pulsar is in a binary with a 7.8-day 
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orbital period with a white dwarf companion, and is thus a 
potential l aboratory for testing the Strong Equivalence Prin¬ 


ciple le.g. iFreire et, all 120121 : lAntoniadis et all 120131 ). The 
integrated pulsed profile of the pulsar consists of a double- 
peaked structure at L-band frequencies with a highly lin¬ 
early polarised trailing component. The amplitude ratio of 
the two components has b een shown to change significantl y 
as a function of frequency dRamachandran fe Kraniei]|2003lf . 
There is also evidence that the ratio is unstable across time 
and may even evolv e on short tim escales of tens of minutes 
dKramer et al.lfl999l ; [P urve 3 Unci). Meanwhile, as the trail¬ 
ing component of the pulse profile is highly polarised, profile 
instability could also arise from improper pol arisation cal¬ 
ibrat ion due to an imperfect receiver model dHotan et al.l 
120041 1 . A better understanding of the profile variation in 
this system is therefore essential in order to further improve 
the timing precision of PSR J1022+1001. This may be used 
to assess and perhaps correct for profile variations in other 
MSPs. 

Examining single pulses has not been common in 
MSPs owing to their typically lower flux densities com¬ 
pared to normal pulsars. Exceptions are the so-called 
giant pulses whi ch have been detected only in a few 
brigh t MSPs (e.g. ICogn ard et all 19961 : iRomani fe Johnstonl 
l200ll : Ijoshi et aid 20041 : Knight et al.l 200 5lf . Analysis of 
the weak single pulses h as been on l y ca rr ied out in 
the f ew brightest MSPs |j enet et ahl 1998^_ Jenet et aid 
200ll: lEdwards fe Stapper3l2003al : Shannon fe Cordesl 20121 : 


Oslowsld_et_al..^014). Based on an autocorrelation analysis, 


Edwards fe Stappersl d2003al l have shown evidence for pulse 


intensity modulation at the phase of the trailing component 
in PSR J1022+1001. However, no detailed investigation into 
the single pulses has been undertaken so far, due to either 
limited system sensitivity or the lack of high-resolution in¬ 
strumentation. 

The rest of this paper is structured as follows: In Sec¬ 
tion [2] we describe the details of observations and data pro¬ 
cessing. Results of the detected single pulses, integrated pro¬ 
file stability, and timing analysis are presented in Section [31 
We conclude and discuss the results in Section [I] 


2 OBSERVATION 

Simultaneous observations of 7 to 9 hours in duration of 
PSR J1022+1001 were conducted with the Westerbork Syn¬ 
thesis Radio Telescope (WSRT) and the Effelsberg 100-m 
Radio Telescope, at four epochs. A summary of the obser¬ 
vations can be found in Tabic [T] The two telescopes were 
chosen as they provide similar sensitivities at L-band so 
that profile instability could be verified from both sites at 
the same time. In addition, Effelsberg is an alt-azimuth 
telescope, meaning that it tracks a source across the sky, 
the changing parallactic angle (PA) results in varied Stokes 
parameters which need to be properly calibrated. On the 
other hand, the WSRT is equatorially mounted with no 
PA variation and thus in principle not influenced by this 
effect. Therefore, using data from these two telescopes al¬ 
lows us to cross-check the result of polarisation calibration 
which was claimed to be the main cau se of profile instability 
dHotan et al.ll2Qo3 : Ivan Stratenll2013h . 

At the WSRT we used the PuMa-II system to record 


Table 1. Parameters of the simultaneous observations of 
PSR J1022+1001 with the WSRT and the Effelsberg 100-m Radio 
Telescope. The symbols T, / 0 , and At s represent the duration of 
observation, overlapping observing frequency range, and the sin¬ 
gle pulse time resolution at the WSRT. 


MJD 

T (hr) 

fo (MHz) 

Af s (ps) 

55909 

7.0 

1300-1440 

4.0 

55974* 

8.4 

1300-1440 

4.0 

56158 

8.5 

1300-1440 

2.0 

56257 

9.0 

1300-1365* 

2.0 


*At this epoch the Effelsberg data were corrupted due to a 
problem with the signal attenuators (see Appendix [A] for more 
details). 

*At this epoch three 25 MHz sub-bands of Effelsberg data were 
lost due to hard drive failure. 

8-bit baseband data dKaruppusamv et alj|200sl j. which were 
later processed offline in two stages. In the first stages the 
8x20 MHz bands were combined to form a multi-band data 
stream corresponding to a 160 MHz contiguous band cen¬ 
tered at 1380 MHz, at the original ti me resolution of 25 ns. 
The DSPSR software (for details see Ivan Straten fe Bailesl 
1201 lT) was used to coherently dedisperse the data over the 
entire bandwidth using a 64-channel synthetic software in¬ 
terbank at the DM of 10.246 cm~ 3 pc and folded to form 
10-s averages in the second stage. Additionally, we searched 
for significant single pulses at this stage. This was done af¬ 
ter summing up the total intensities of all frequencies, and 
averaging by the number of frequency channels. Any pulse 
period that contained a peak greater than five times the 
root-mean-square (rms) of the off-pulse noise was written 
out as a single pulse candidate. In the post-processing stage, 
these were then confirmed to be true pulses by comparing 
the phase to that of the integrated pulse profile. 

At Effelsberg, the new PSRIX pulsar backend (Karup- 
pusamy et al., in prep) was used in baseband mode. These 
data were reduced in the same way as the WSRT data, ex¬ 
cept that here, 8x25 MHz bands were produced with a cen¬ 
tral frequency of 1347.5 MHz. In addition, at regular inter¬ 
vals of 40 min, the telescope was pointed 0.5° off-source for 
90 s, when a pulsed noise signal was injected at 45° into the 
feed probes as a calibrator for polarimetry. 

The 10-s integrations and single pulses w ere then pro¬ 
cesse d with the psrchive software package dHotan et al.l 
l2004ll . For the 10-s integrations from Effelsberg we per¬ 
formed polarisation calibration with the single-axis receiver 
model which corrects for the differential gains between the 
two feedfl Next we used pazi (psrchive’s interactive 
RFI zapper) to clean the radio frequency interference (RFI) 
by visual examination of the averaged pulse profiles and sin¬ 
gle pulse candidates. Finally, we selected all 10-s integra¬ 
tions that overlapped in time between the two telescopes, 
and then removed the non-overlapping frequency channels, 
leading to an overlapping observing frequency range of 
1300 - 1440 MHz. 

Example polarisation profiles of PSR J1022+1001 
from both the WSRT and the Effelsberg observation at 
MJD 56257 centred at 1320 MHz, can be found in Fig. |T] 


1 http://psrchive.sourceforge.net/manuals/pac/ 
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Figure 1 . Polarisation profiles of PSR J1022+1001 obtained at 
MJD 56257 from the Effelsberg (middle) and the WSRT (bottom) 
observation at 1.3 GHz. The P.A. in the top panel stands for the 
polarisation position angle of the linear component. The Effels¬ 
berg and the WSRT data show consistent P.A. curves and here 
we plot the one from the Effelsberg data as a demonstration. The 
solid, dashed and dotted lines represent the Stokes parameters I, 
L and V, respectively. There is a “negative” dip, i.e., an under¬ 
estimation of system temperature, on the WSRT total intensity 
profile (at phase 0.5 2) which corresponds to t he low-bit sampling 
artefact discussed in ljenet fc Andersonl (119981) . This is due to the 
fact that the signal was originally quantized using 2 bits per sam¬ 
ple at each individual dish of the WSRT. Further investigation 
shows that the dip is mostly from high frequency channels within 
each 20 MHz band. The total intensity profile is slightly lowered 
by this underestimation, which makes the valley region between 
the two components appear over polarised. Further investigations 
can be found in Appendix [A] 


The integrated profile shows a double-peaked structure, with 
the trailing component consisting of highly linearly polarised 
emission. The polarisation properties are in agreement be¬ 
tween the tw o sites, with a c orrelation coefficient of 0.999 
(as defined in iLiu et al ] E5a . Given the profile instability 
of this pul sar, especiall y in observations separated by sev¬ 
eral years (|Purverll201(t ) . the polarisation com ponents quali¬ 
tatively agree with the previous observations l|K ramer e t al.l 
1 19991 : Ramachandran fc Kranierll2003l : Ivan Stratenll201.il '). 


3 RESULTS 

3.1 Sub-pulse properties 

In a total of 35 hours of observations, we have approximately 
14,400 single pulse detections above our 5-<r threshold, which 
corresponds to about 700 sub-pulsefl located at the lead¬ 
ing component of the integrated profile and about 13,700 at 
the trailing component. The detections were all made from 
the observations at MJD 55909 and 56257 as the pulsar 


2 In each rotation, the pulsar produces a single pulse which can 
be composed of one or more sub-pulses. This better reflects our 
situation as occasionally, more than one sub-pulses can be de¬ 
tected from a single period. 

3 The leading and trailing component pulses are distinguished by 
an arbitrary bound of rotational phase 0.485 in Fig. |T| 


Figure 2. Sub-pulses detected at MJD 56257 as a function of 
time and their peak signal-to-noise ratios (S/N s ) divided by the 
local average S/N of all single pulses observed during the 5-min 
interval centred on the epoch of the single pulse ((S/N s )5 m i n ). 
The location of peak bin was selected based on a 1 6-//S time reso¬ 
lution. The top panel shows detections at the phase of the leading 
component in the averaged profile, and the bottom shows those 
at the phase of the trailing one. 


was significantly weaker due to interstellar scintillation at 
MJD 55974 and 56158. The highest peak signal-to-noise ra¬ 
tio of detection is close to 10. Fig. [2]shows the sub-pulse de¬ 
tections achieved at MJD 56257 at their occurrence epochs 
and peak signal-to-noise ratios (S/N s ) relative to the local 
average S/N of all single pulses observed during the 5-min 
interval centred at the epoch of the single pulse ((S/N s )5 m i n ). 
We can see that there is an evolution in the ratio which is 
due to the pulsar flux decreasing because of interstellar scin¬ 
tillation. With our sensitivity, useful detections are feasible 
only when integrating 5-min observations results in a profile 
of peak S/N (S/N5 m i n ) above 30. The detection rate is high¬ 
est between MJD 56257.055 and 56257.060, corresponding 
to 4.5% of all rotations. We have calculated time separations 
between neighbouring sub-pulses at the trailing component 
within this period of time, and over 95% of them are less 
than 1 s, indicating that pulse nulling with duration longer 
than Is (~ 60 rotations) is not detected. 

Fig. [3] shows the distribution of the peak amplitude 
phases for the detected sub-pulses, compared with the in¬ 
tegrated profile. All detections were obtained within the on- 
pulse phase and most of them are clustered within the range 
defined by the trailing component in the average profile. 
The two maxima correspond to the two peaks in th e phase- 
resolved modulation index of lEdwards fc Stappera (|2003bT ). 
In Fig. [4] we plot the distribution of the flux densities of 
the sub-pulses divided by the averaged flux densities of a 
single pulse in 5-min integrations. The plot demonstrates 
that none of the detections has a flux density that is above 
3.8 times the average, and there is therefore no evidence for 
giant pulses. 

Fig. [5] presents example polarisation profiles of sub¬ 
pulses at the phase of the leading and trailing components. 
It can be seen that pulses from the leading component do 
not show significant polarisation, while those from the trail¬ 
ing component are highly linearly polarised. The polarisa- 
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Figure 3. Distribution of peak amplitude phase for all detected 
sub-pulses (lower panel), in comparison with the integrated pro¬ 
file (upper panel). Again, the peak bins were chosen based on a 
16-^s time resolution. The integrated profile is from the WSRT 
observation at MJD 56257. 



Pulse flux / mean flux 


Figure 4. Distributions of relative flux densities to averaged sin¬ 
gle pulse flux densities in 5-min integrations, for sub-pulses de¬ 
tected at the leading and trailing components, respectively. No 
pulse with flux above 3.8 times the average has been seen. 


tion fractions of both components are consistent with the 
integrated pulse profile, suggesting that the detected sub¬ 
pulses represent the bright end of the energy distribution 
of all pulses. The bottom plot in Fig. [5] shows the averages 
of the leading and trailing component pulsefl respectively, 
whose polarisation characteristics are in agreement with the 
corresponding individual detections. It is interesting to note 
that in both averages, the components that do not include 
detected sub-pulses are consistent with the average of all 
periods in total intensity. This means that during the occur¬ 
rence of bright emission in one component, there still exists 


4 Here we selected pulses with 5 —a detection in only one compo¬ 
nent. Occasionally, bright emission can occur at both components 
within a single period, which will be discussed in the next para¬ 
graph. 


emission at an average level at the location of the other 
component. 

Meanwhile, sub-pulses can be detected from the leading 
and trailing components within the same pulse period. In to¬ 
tal, about 280 such events have been detected. In all events, 
the two phases of peak amplitude are well separated from the 
phase bound 0.485 (used to group the detections), exclud¬ 
ing detections of a single pulse with broad width spreading 
over the region of both components. We will later refer to 
such an event as a ‘bi-pulse’, and their corresponding inte¬ 
grated profile is shown in Fig. [S] The valley region between 
the two components is significantly steeper in this profile 
than the ordinary average in Fig. [l] while the polarisation 
properties of the two components are consistent. To investi¬ 
gate the correlation of occurrence between the leading and 
trailing component pulses, we carried out a further statis¬ 
tical study. Pulses were selected at both epochs, from the 
period of time when the leading component pulses could be 
detected. A summary of the selections can be found in Ta¬ 
ble^ We formed the statistical test problem in the following 
statement: 

{ Ho Leading and trailing pulses are uncorrelated. 

Hi Leading and trailing pulses are correlated. 

The bi-pulse phenomenon is defined as: both the leading 
and trailing component pulses are detected within a single 
period. Thus, under the null-hypothesis Ho, in each pulse 
period the probability of observing bi-pulses is 

Pb = PPt , (1) 

where P\ t , P\ and P t are the probabilities of observing bi, 
leading, and trailing component pulses, respectively. By 
combining the statistics from these two epochs, we have 
calculated the observed Pi and Pt which gives an expected 
value of Po as 3.80 ±0.16 x 10~ 5 under the null-hypothesis 
Ho- Meanwhile, the observed bi-pulse occurrence probabil¬ 
ity was calculated to be 6.32 ± 0.39 x 10~ 4 , about 17 times 
larger. For a more quantitative study, over N pulse periods 
the statistical distribution function / of the total number of 
bi-pulses, leading and trailing component pulses are 


f(N h \P h , N ) = C% b P? b ( 1 - P h ) N ~ Nh , 

(2) 

f(Ni\P u N) = Cglf'i 1 - Pi ) N ~ N ', 

(3) 

f(N t \Pt,N) = C^P t Nt ( 1 - Pt) N ~ Nt , 

(4) 


where C% is the binomial coefficient defined as CJ) = 
(fe!)/(a!(6—a)!). Using the Bayesian theorem with a flat prior 
for the distributions of probability of each type of pulses, we 
have 

f(Pi\N h N) = f(Ni\P u N)(N + l), (5) 

f{Pt\N t ,N) = f(N t \Pt,N)(N + 1), (6) 

where f(Pi\Ni,N) and f(Pt\Nt,N) are the posterior for Pi 
and Pt. Thus the probability distribution of Nh given Ni 
and Nt is 

f{N b \Ni,N t ,N) = [ [ dPidPtC% h (PiP t ) Nh 

Jo Jo 

(1 - plp t ) N - Nb 
f(Pi\N h N)f(Pt\N t ,N). 


( 7 ) 
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Figure 5. Example polarisation profiles of sub-pulses detected 
at MJD 56257, at the phase of the leading (top) and trailing 
(middle) component in the average profile. The bottom plot shows 
the integrations of detected sub-pulses from the leading (upper 
panel) and trailing (lower panel) components, respectively. The 
solid, dashed and dotted lines stand for the Stokes I, L and V, 
respectively. The dashed-dotted line represents the averaged total 
intensity of all periods. 


Figure 6. Polarisation profile averaged over all detected bi¬ 
pulses, where sub-pulses were detected at both components within 
a single period. 

Table 2. Summary of selected pulses used to investigate the oc¬ 
currence correlation between the leading and trailing component 
pulses. Here N, Np Nt and Nb denote the number of periods 
within the time range the pulses were chosen from, number of 
leading and trailing component pulses, and number of bi-pulses, 
respectively. The probability of detecting more than TVb bi-pulses 
is denoted by P(n ^ TV i) 


Phase 

MJD 55909 

MJD 56257 

N 

2.1 x 10 5 

2.0 x 10 5 

Ni 

72 

525 

N t 

2737 

7929 

N b 

25 

234 

P(n >= N b )* 

io - 252 

10-88! 


* Due to the difficulties in the numerical evaluation of the 
integral, we used asymptotic values for the error function. This 
could result in one order of magnitude error, but will not change 
the conclusion. 

Under the null hypothesis Ho, we have 

N B 

P(N h >N B orN h <N A ) = l~Y / f(N h \N u N t ,N), (8) 

n a 

where P(Nb > Nb or JVb < Na) is the probability of find¬ 
ing that the occurrence of bi-pulses Nb is greater than the 
threshold IVb or less than the threshold N A , under the as¬ 
sumption that the leading and trailing component pulses 
are uncorrelated. This provides a good statistic for testing 
whether the observations agree with the null hypothesis Ho. 
Accordingly, we have estimated the probability to detect 
more than Nb bi-pulses (P(n >= Nb)) based on the statis¬ 
tics in Table[2] which appears to be very close to zero at both 
epochs. We therefore conclude that the leading and trailing 
component pulses are correlated in their occurrence, other¬ 
wise it would be extremely rare to observe such a high value 
of Nb with the null-hypothesis Ho . In young pulsars, corre¬ 
lated emission modulation at different rotational ph ases has 
been found in a few cases llWeltevrede et ahll2012l . Karup- 
pus amy et al. in prep.) . 

Ijenet et al.l (I1998T) showed an inverse correlation be¬ 
tween pulse width and peak amplitude for the sub-pulses of 
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Figure 7. Distribution of pulses from the trailing component 
by their relative S/N s and width. To avoid bias caused by the 
time-variant detection threshold, only pulses with epoch between 
MJD 56257.042 and 56257.066, and (S/N p )/((S/N a ) 5min ) > 7, 
were selected. Again, the peak bins were selected based on a 16- 
ps time resolution. 

PSR J0437—4715. In Fig. □ we perform the same investiga¬ 
tion by showing the distribution of detected pulses from the 
trailing component on a width—relative S/N s plane. Here, to 
avoid bias caused by a variable detection threshold, we only 
selected detections between MJD 56257.042 and 56257.066 
and with (S/N s )/((S/N s )5 m i n ) > 7. It is interesting to note 
that, instead of an anti-correlation, there exists a favoured 
pulse widtl0 value of about 0.25 ms which is consistent for 
pulses of different relative S/N s values. The probability of 
occurrence also does not favour pulses with narrow width 
and high peak amplitude. Fig. [8] gives the distribution of 
the same pulses with respect to their phases and relative 
S/N s . The plot demonstrates that the occurrence phases of 
the sub-pulses are strongly clustered, especially for those of 
peak 10 times above average. The phenomenon that brighter 
p ulses show up earlier i n pha se, found for PSR J1713+0747 
in I Shannon fc Cordesl J20121 ). is not seen. The property of 
regular shape and phase could possibly be translated into 
precise pulse time-of-arrival (TOA) measurements, which 
will be further investigated in Section U.31 

3.2 Variation of integrated profiles and sub-pulse 
occurrence 

As mentioned in Section |T| it has not been conclusively 
shown whether the average profiles of PSR. J1022+1001 ex¬ 
hibit intrinsic shape variation. Here, fo llowing the idea of 
previous analyses (e.g. iKramer et al.lll999l ) , we calculate the 
amplitude ratio of the leading and trailing components of to¬ 
tal intensity profiles for 10-min integrations, and plot them 
against time. For each component, the amplitudes were de¬ 
fined as the total intensity within a fixed phase range (1.7% 
of a period), after subtracting the baseline. Their errors were 
calculated from the rms of the off-pulse bins. We used a 

5 Here the pulse width is defined by the on-pulse region. 


Figure 8. Distribution of the same pulses as in Fig. [71 on a 
phase—relative S/N s plane. 

40 MHz sub-band ranging from 1300 to 1340 MHz, where for 
most of our observing time, the flux density is the highest 
within the entire overlapping band. The results from data 
at MJD 55909 and 56257 are shown in Fig. [9] where at 
both epochs the measurements show a clear evolution with 
time and are consistent between the two sites. The variation 
tren d is simil ar to what has been indicated in Fig. 3.6 of 
iPurved (Mh . We also calculated the weighted correlation 
coefficients ( p w s) between trends from the two sites, which 
is defined as 

COV w (£i, Wi, crf^i, <j w i ) 

Pw — , • 

^/covw (£ i , £i, as , i , ae , i ) cov w (W;, W;, a w ,;, aw , ;) 

(9) 

Here £i, W; are the ratio measurements from Effelsberg and 
WSRT, respectively, a w,i are their measurement er¬ 

rors, and cov w is the weighted covariance defined by 

^2(xi - x){yi - y)/(o’x,i a y, i ) 

COV w (r,, yi, &x,ii &y,i) — ^ ^ , 

( 10 ) 

where the overline stands for weighted mean. The calcu¬ 
lated p w is 0.60 at MJD 55909 and 0.91 at MJD 56257, 
respectively, confirming that measurements from the two 
sites exhibit the same evolution trend. Note that at the be¬ 
ginning of MJD 56257 the pulsar was significantly brighter 
and the profile was showing a steeper variation trend than 
at MJD 55909, which thus corresponds to a value of p w 
closer to unity. To understand the estimation error of p w , 
we performed a Monte Carlo simulation. In each iteration, 
we randomly altered the measured component amplitudes 
based on Gaussian distributions with the standard devia¬ 
tion equal to the measurement errors, and thus calculated 
a new amplitude ratio and p w - Fig. HOI shows the histogram 
of p w with 10,000 iterations obtained from measurements 
at MJD 56257, where p w is larger than 0.78 and 0.64, with 
68% and 95% confidence level, respectively. The same inves¬ 
tigation has also been performed for the measurements at 
MJD 55909 and the corresponding thresholds are 0.40 and 
0.20, respectively. These mostly rule out the possibility of a 
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MJD-56257 

Figure 9. Amplitude ratio measurements of the two compo¬ 
nents in the PSR J1022+1001 total intensity profiles from two 
epochs. The measurements were made using 10-min integrations 
of a 40 MHz band centred at 1320 MHz. Note that measurements 
of fractional error larger than 10% were discarded, which leads to 
the gap within the results at MJD 56257. 


Table 3. Weighted correlation coefficient (p w ) between the am¬ 
plitude ratio measurements from the WSRT and Effelsberg data 
at two epochs, as well as their lower limits with 68% (£68) and 
95% (£ 95 ) confidence level. 


MJD 

Pw 

£68 

£95 

55909 

0.60 

0.40 

0.20 

56257 

0.91 

0.78 

0.64 


non-correlation between the two sites. A brief summary of 
the results can be found in Table [3 

To investigate the effect of polarisation calibration on 
our results, we also applied the template-matching calibra¬ 
tion te chnique to the E ffelsberg data at MJD 55909 and 
56257 (jvan Stratenll200~6l . Lee et al. in prep.), based on the 
template polarisation profi le from Europea n Pulsar Network 
Pulse Profile Database (IStairs et al.lll999T) showing consis¬ 
tent polarimetry with previous published results. For this 
calibration scheme, we used both the single-axis and the re¬ 
ception description of the receiver properties, the latter of 
which models bot h the differentia l gain and cross-co upling 
of the two feeds llOid et al.ll2004l ; Ivan Stratenl l2003 l. The 



Figure 10. Distribution of 10 4 simulated values of p w based on 
amplitude ratio measurements from data obtained at MJD 56257. 


calibrated data were then processed by following the same 
procedure as described above. Both receiver models led to 
a fit with reduced \ 2 close to unity. Still, the effect of feed 
cross-coupling was estimated to influence the total inten¬ 
sity profile by less than 1%, and the true value could not 
be measured due to limited sensitivity. As a demonstration, 
in Fig. [TT]we show the obtained polarisation profiles from 
MJD 56257 data with the calibration result based on cali¬ 
brators and the single-axis model shown in Fig. [T] It can 
be seen that the three profiles from different calibration 
schemes or models exhibit consistent polarisation proper¬ 
ties, especially the two achieved with the single-axis model. 
Using the reception model results in a slightly lower linear 
polarisation percentage (middle panel), mostly due to ab¬ 
sorption of power into the feed leakage terms. In addition, 
the calculated p w s between the new amplitude ratio mea¬ 
surement trend and the WSRT result, is 0.89 when using 
the reception model and 0.91 when single-axis calibration 
was applied. These are almost the same as the value ob¬ 
tained above. The same analysis at MJD 55909 data led to 
p w = 0.55 when using the reception model and p w = 0.60 
when applying single-axis calibration, also consistent with 
the result based on calibrators. Further investigation based 
on independent Effelsber g data for the Large Europ ean Ar¬ 
ray for Pulsars project (IKramer fe Stappersl 1201011 . shows 
that the circular receiver does not suffer from strong cross¬ 
coupling between the two feeds (Lee et al. in prep.), which 
validates the use of the single-axis model in this case. There¬ 
fore, it is highly unlikely that our detected profile variation 
is due in large part to improper polarisation calibration. 

Note that the profile shape has a clear fr equency de¬ 
pendence (e.g. iRamachandran fe Krarnerl 120031 ). If the flux 
density distribution within the observing band changes sig¬ 
nificantly in time due to ISM diffractive scintillation, the av¬ 
eraged profiles (after frequency summa tion) may there fore 
also exhibit time-variant behaviour (e.g. iLiu et alll201l| j. To 
understand its contribution to our results, here we simulate 
the varying trends by this effect and compare them with 
the real measurements. At each epoch, we first used an in¬ 
tegrated profile for the entire duration of the observation to 
establish a profile model with frequency dependent shape. 
As an example, Fig. [12] shows the measured amplitude ratios 
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Figure 11. Comparison of Effelsberg polarisation profiles ob¬ 
tained from data at MJD 56257, based on calibrator+single- 
axis model (left), template+reception model (middle), and 
template+single-axis model (right). The solid, dashed and dot¬ 
ted lines stand for the Stokes parameter f. L and V, respectively. 
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Figure 12. Measured amplitude ratios as a function of frequency, 
from the profile models obtained at MJD 55909. The measure¬ 
ments from the Effelsberg and WSRT data are based on profiles of 
bandwidth 6.25 and 5 MHz, respectively. The edges of the 20 MHz 
sub-bands (starting from 1300 MHz) in the WSRT data have been 
marked by dashed lines. Note that for the WSRT data there is an 
additional increase of the measurement values within a sub-band, 
associated with the low-bit digitisation artefact shown in Fig. [T] 
More details of this effect are shown in Appendix [Al 


as a function of frequency^ from the profile models for the 
MJD 55909 data. We then applied the measured frequency- 
dependent flux densities from each 10-min integration to the 
model to create a simulated profile, and computed the re¬ 
sulting amplitude ratio after frequency summation. For this 
purpose, we kept the frequency resolution of 6.25 MHz for 
the Effelsberg data and 10 MHz for the WSRT data, so as 
to have enough signal within each frequency channel. The 


6 The defined phase ranges to calculate component amplitudes 
are identical for all central frequencies as the data were already 
de-dispersed as mentioned in Section [2] 


channel width is significantly smalle r than the estimated 
scintillation ban dwidth of this pulsar dCordes fe Lazioll2002l : 

I You et aLll2007l) . The uncertainties in the amplitude ratios of 
the simulated profiles were calculated from the measurement 
errors of the flux densities estimated from each frequency 
channel of each 10-min integration. The flux densities as 
well as their errors were measured by using pdv (psrehive’s 
archive data displayer), where the flux density is defined as 
the total intensity within a pulse period after subtracting 
the baseline, and their errors were calculated from the rms 
of off-pulse phase bins defined by 3-cr threshold. By follow¬ 
ing this procedure, the observed variation trend would be 
reproduced if the profile instability is mostly dominated by 
the scintillation effect. In Fig. [13] we present the results at 
MJD 55909 and 56257 from both sites, based on both the se¬ 
lected 40 MHz sub-band and the entire bandwidth. Clearly, 
in most cases the simulation resulted in a flat trend and did 
not reproduce the actual measurement^. The reduced x 2 s of 
the measurement series with respects to their corresponding 
simulated trenejf] fell in between 2.0 and 22, all significantly 
above unity. We therefore conclude that our detected profile 
instability shown in Fig.[9]is not significantly affected by the 
scintillation effect. 

We do not detect profile variation at the other two 
epochs, MJD 55974 and 56158, mostly because the source 
was comparatively dim due to scintillation. At these two 
epochs, with the WSRT 160 MHz observations the S/Ns m i n 
values are all below 30, while the maxima at MJD 55909 
and 56257 reach 75 and 110, respectively. This may explain 
the contradictory conclusions from previous studies on pro¬ 
file stability on short timescales. In addition, it can be seen 
from Fig. [9] that the amplitude ratio variation may not be 
significant enough for detection within 0.1 day, which is the 
longest observation used in previous work. Therefore, detec¬ 
tion of profile instability is achievable only when the pulsar 
is bright and the variation is large. 

We note that our analysis draws seemingly different 
conclusi ons on the profile stability of PSR J1022+1001 com¬ 
pared to lLiu et al.l d2012l) . due to a few reasons. Firstly, our 
new observations detected the pulsar when it was signifi¬ 
cantly brighter than in the previous ones. For instance, at 
the beginning of the MJD 56257 observation an integration 
time of 1800 s led to a profile with signal-to-noise ratio over 
800, while integrat in g the e ntire observation (approximately 
1,900 s) in I Liu et al.l d2012l) only result e d in a value of about 
480. Secondly, the work in iLiu et all (l2012h was based on 
1-min integrations and thus focused more on searching for 
random variations in profile shapes on short timescales. Fi¬ 
nally, the duration of the previous observation does not pro- 


7 Note that an exceptional case may be drawn from the simu¬ 
lation on full bandwidth of the WSRT data from MJD 55909. 
This indicates that profile evolution based on large bandwidth 
could be significantly contributed by scintillation effect. The rea¬ 
son why the simulation based on Effelsberg data from the same 
epoch did not reproduce the variation trend, is that our Effels¬ 
berg profile shows significantly less frequency dependency and is 
thus less affected by the scintillation effect. 

(ri - Si ) 2 


8 Defined as 


N r — 1 . 


r z . (j z ■ 
r,i ' s,i 


where N r is the number of 


amplitude ratios, r^, Si are amplitude ratios from real measure¬ 
ment and simulation, respectively, and cr r ,i , &s,i are their errors. 
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Figure 13. Simulated peak amplitude ratios based on flux densities measured from 10-min integrations, compared with the actual ratio 
measurements. The simulations were carried out on both the WSRT (left) and Effelsberg (right) data at MJD 55909 (top) and 56257 
(bottom), using both the 40 MHz sub-band and the entire available bandwidth. The uncertainties of the simulated ratios are typically of 
order < 10 3 , and thus not visible from plots. Note that for MJD 56257, the analysis was carried out only for the first half of the data, 
as in the second half the flux densities within individual frequency channels were not measurable on many occasions due to low S/Ns. 
The effective observing bandwidth of Effelsberg at MJD 55909 is 160 MHz due to a receiver cutoff below 1280 MHz. 


vide the sensitivity to profile variations on a timescale of 
several tens of minutes detected in this paper. 

In principle, the profile variation can be fully attributed 
to the detected bright sub-pulses, if either their shape or oc¬ 
currence rate varies by a sufficiently large amount in time. 
Such a possibility has been investigated in Fig. 1141 To bal¬ 
ance between potential shape bias introduced by a vari¬ 
able detection threshold and time span allowing us to see 
profile variation, here we selected sub-pulses detected with 
(S/N a )/((S/N s ) 5m in) > 9 and MJD between 56257.042 and 
56257.074. It can be seen that while the amplitude ratios 
from 5-min integrations is gradually increased in time (« 6% 
from MJD 56257.04 to 56257.07), no similar variation trend 
is witnessed from the averages of the sub-pulses. This in¬ 
dicates that the profile variation is not caused by the se¬ 
lected bright sub-pulses. In fact, the number of detected 
sub-pulses is approximately 200 — 300 per 5-min window 
(1 — 1.5% of all periods). Therefore, the amplitude ratios 
of the sub-pulse averages would need to change by a fac¬ 
tor of 4 or more if variability in the 1 — 1.5% of detected 
sub-pulses were the cause of the profile variation in 5-min 
integrations. Note that the number of detected sub-pulses 
remain consistent in time, though it is proportional to the 


S/N of 5-min integrations, which is expected due to vari¬ 
ability in detection sensitivity. This suggests that the profile 
variation is not induced by the occurrence rate variation of 
bright sub-pulses. In fact, assuming profile variation mainly 
occurs at the trailing component, the increase of amplitude 
ratio from 5-min integrations corresponds to 6% decrease 
of the flux density at the trailing component. Given the 
average amplitude ratio of 0.18 observed in the detected 
sub-pulses and an amplitude ratio of 1 for ordinary inte¬ 
grations, the selected sub-pulses contribute 6-9% of the flux 
density at the trailing component. Therefore, if the change 
of occurrence rate of the selected sub-pulses were fully re¬ 
sponsible for the amplitude ratio variation, the number of 
detections around MJD 56257.07 would have dropped by 
70-100%, which is clearly not the case. Therefore, it is un¬ 
likely that our selected sub-pulses cause the detected profile 
variation, which suggests that sub-pulse variability and in¬ 
tegrated profile variation are two different effects. Neverthe¬ 
less, the analysis is still restricted to roughly the brightest 
1% of all single pulses, and will be significantly improved by 
increasing system sensitivity. 
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Figure 14. Top to bottom: amplitude ratio measurements from 
5-min integrations and from sub-pulse averages, numbers of sub¬ 
pulses in 5-min, S/Ns of 5-min integrations. Here the sub-pulses 
were selected with (S/N B )/((S/N s ) 5 m ; n ) > 9 and MJD between 
56257.042 and 56257.074. The 5-min integrations were formed 
from WSRT data. 

3.3 Timing of sub-pulses 

To investigate the impact of sub-pulses on pulsar timing, 
we used the sub-pulses selected in Fig. [7] (3.4% of all pe¬ 
riods, from the trailing component, in total approximately 
4,500 pulses) and compared their timing with ordinary inte¬ 
grations based on all rotations. The TOAs are measured by 
foll owing the cl assic template-matching method described 
in iTavlod Il992h . The template profile for timing based on 
sub-pulses is formed by summing all detected sub-pulses 
and perform i ng a Gaussian component smoothing as in 
e.g. iKrameri (Il994h . The timing resid uals were calculate d 
with the TEMP02 software package dHobbs et al.l l2006h . 
Here we used the timing solution obtained from the Euro¬ 
pean Pulsar Timing Array collaboration (Desvignes et al. in 
prep.), without fitting for any parameters. 

In Fig. [15] we present the timing residuals calculated 
with averages of 30 pulses (not contiguous in time), equiva¬ 
lent to an integration time of 0.5 s. The timing solution has 
an rms of 4.3 /xs which is greater than the average of the 
errors owing to white noise (1.6/xs). This suggests the exis¬ 
tence of phase jitter. To better understand the timing noise, 
a standard Kolmogorov-Smirno v (K-S) t e st ha s been per¬ 
formed on the TOAs as done in iLiu et ah! ll2012h . The mea¬ 
sured p-value (~ 0.98) is close to 1, indicating no significant 
deviation of the r esiduals from a Gaussian distribution (e.g. 
[Press et ahlll986l 'l. In Fig. [16] we divide the entire band into 
two sub-ban ds, and present a corre latio n plot of TOAs from 
them, as in lOslowski et al] d201 lh and IShannon fc Cordesl 
ll2012f) . The TOA pairs are shown to be highly correlated 
and the rms of each individual group is also estimated to be 
roughly equal (4.6/xs), meaning that the noise is correlated 
and dominated by phase jitter. 

In contrast, we did not detect any significant impact of 
pulse phase jitter on timing with ordinary integrations. Data 
obtained between MJD 56257.044 and 56257.080, when the 
observing S/N was maximal, have been examined based on 
1-min averages. The resulting timing rms is 1.3/xs, close to 
the estimated white noise level (1.1 /xs) from the TOA errors. 


Figure 15. Timing residuals based on 151 averages of 30 sub¬ 
pulses. The rms is 4.3 /xs while the white noise rms indicated from 
the TOA errors is 1.6 /xs. 



At, (ps) 

Figure 16. Correlation plot of TOA pairs calculated based on 
sub-pulse averages from two individual 80 MHz sub-bands. The 
data used here are the same as in Fig. [El] 

Subtraction of the white noise contribution from the timing 
residuals leads to an upper limit of ~ 700 ns for jitter noise 
based on a 1-min integration time. 

In Fig. [lT]we present the timing rms values when alter¬ 
ing the number of integrated sub-pulses and compare them 
with timing rms yielded by ordinary integrations. It is shown 
that the jitt er noise decreases followi ng the crj oc tint law 
as expected dCordes fe Shannonll201(ih . where tint is the in¬ 
tegration length of each average. Averaging 30 sub-pulses 
leads to 151 TOAs of rms 4.3/xs. Considering that during 
the same period of time ordinary 1-min integrations lead to 
34 TOAs of rms 1.3/xs, timing the averages of sub-pulses 
is almost as effective as timing the ordinary integrations 
(4.3/^151/34 ~ 2.0/is, close to 1.3/xs). Note that only 3.4% 
of the single rotations were picked up, an improved system 
sensitivity would enable significantly more sub-pulse detec¬ 
tions. Including those sub-pulses in timing may significantly 
decrease the jitter noise and make the timing based on sub¬ 
pulses more efficient (achieving more measurements of the 
same rms in the same observing time, or equal number of 
































































































































Profile instability of PSR J1022+1001 11 


10 


C/3 

C/3 

E 

03 

O 


1 

0.1 1 10 100 
Integration time (s) 

Figure 17. Timing rms of averages based on different numbers of 
bright sub-pulses, compared with rms from integrations of all pe¬ 
riods. The dashed line shows a power-law modelling with a fitted 
index of — 0.49T0.02. Note that an integration of 10s corresponds 
to roughly 600 pulses. 
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measurements with lower rms) than using all-period integra¬ 
tions. However, this extrapolation is based on the assump¬ 
tion that the sub-pulses to be included exhibit the same 
phase jittering as the currently selected ones, which requires 
further investigations with improved sensitivity to validate. 


4 CONCLUSION AND DISCUSSION 

In this paper, we have performed a detailed investigation 
into the profile variation of PSR J1022+1001. Within our 
35-hr observations, approximately 14,400 sub-pulses were 
detected above a 5 -a threshold from the WSRT baseband 
data, and most of them are coincident with the trailing com¬ 
ponent of the integrated profile. The flux densities and po¬ 
larisation properties suggest that they are the bright end 
of the pulse energy distribution and not a separate pop¬ 
ulation. No giant pulse has been detected. Pulses from the 
leading and trailing components can happen within the same 
rotation and their occurrences are shown to be correlated. 
For pulses from the trailing components, a preferred pulse 
width of approximately 0.25 ms has been found. Using si¬ 
multaneous observations at the Effelsberg telescope and the 
WSRT, we have tracked the variation of integrated profiles 
on a timescale of several tens of minutes. Profile instabil¬ 
ity due to improper polarisation calibration and diffractive 
scintillation was excluded. We have not yet discovered an 
association between integrated profile instability and sub¬ 
pulse occurrence, possibly owing to the present sensitivity 
limit. In addition, we have demonstrated the dominance of 
phase jitter in timing residuals with the sub-pulse averages 
from the trailing component, and obtained an upper limit 
of ~ 700 ns for jitter noise based on continuous 1-min inte¬ 
grations. Further investigation with better sensitivity is still 
required to find out if timing the averages of sub-pulses is 
more efficient than using continuous integrations. 

The spin-down power (E) of PSR J1022+1001 is esti¬ 
mated to be 3.8 x 10 32 ergs s -1 (obtained from AT NF Pul¬ 
sar Catalog, for details see [Manchester et al.l[2005h . signifi¬ 


cantly less than t hose of the current f ew MSPs with giant 
pulse discoveries (|Cognard__elLaL 19961: Romani fe Johnstonl 


l200ll : l.loshi et al.1120041 : iKnight et al.l 20051 ). Thus. our non¬ 
detection of giant pulses from PSR J1022+1001 tends to 
coincide with the correlation between giant pulse emissiv- 
ity a nd spin-down luminosity in MSPs (e.g. IKnight et al.l 
l2005lf . Still, the detected sub-pulses satisfy the criterion of 
the so-called “giant micropulses”, which requires the peak 
flux densitiy much greater than 10 times that of the aver¬ 
age profile, but the i ntegrated flu x density not exceeding 
10 times the average dCairnafeoOlf ) . Our detection, togethe r 
with the discovery in PSR J0437—4715 d.Tenet et al.lll998l ). 
suggests that emissivity of giant micropulses in MSPs is not 
necessarily related to spin-down luminosity. This suggests 
the possibility that giant micropulses may have a different 
origin from giant pulses. 

The two components of PSR J1022+1001 are separated 
by only 5% of the period, but have dramatically different 
emission behaviours. It could be that the two components 
actually originate from well-separated regions in height, du e 
to the strong winding of the magnetosphere fe.g. |Petil]|201l| j. 
A fit of the aberration & retardation field-line model to the 
integrated profile which estim ates the emission latitude, may 
shed more light on this issue (iDvks et al.ll2003 '). 

Admittedly, single pulse studies of MSPs are still mostly 
limited by detection sensitivity. With the next generation of 
radio telescopes, such as the Five-hundred-metre Aperture 
Spherical Radio Telescope and the Square Kilometre Ar¬ 
ray, the system sensitivity for pulsar observations will be 
incre ased b y 1 -2 ord ers of magnitude from the current level 
fe.g. ILiu et al.ll201in . This will enable more comprehensive 
investigations into MSP single pulses, by increasing the num¬ 
ber of detections, improving the detection qualities, and en¬ 
abling more case studies. Such efforts would be greatly help¬ 
ful in understanding the behaviour of single pulse emission 
as well as its impact on pulsar timing stability, and thus de¬ 
veloping potential approaches to improve timing precision. 
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Figure Al. Gray-scale image of the profile against frequency 
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time. 
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APPENDIX A: EXAMPLE OF CORRUPTED 
DATA DUE TO FAULTY POWER LEVEL 
SETTING 

The data from the failed observing session at MJD 55974 
at the Effelsberg telescope provide insights into the effect of 
2-bit sampling on profile evolution. In this observing run, 
the signal level attenuator was inadvertently set too high, 
resulting in voltages that were too small as input to the 8- 
bit analog-to-digital converter (ADC) and reducing it to a 
low-bit (2-bit or even smaller) system. This data allow us 
to assess the effect of low-bit sampling on PSR J1022+1001 
profile shape variation, thus highlighting an extreme case of 
signal clipping. 

Fig. lAll shows a gray-scale image of the profile as a func¬ 
tion of frequency before dedispersion. Within each 25 MHz 
sub-band, there is negative po wer distributed in a simila r 
fashion as shown by Fig. 4 of Ijenet fe Andersonl (1l998lh 
which is a consequence of low-bit sampling. As the pulsar 
signal is dedispersed, the negative dip which would suppress 
the power of the signal, leads the profile at lower frequencies 
and trails at higher ones. This will result in additional peak 
amplitude ratio variation across the sub-band which can be 
witnessed in Fig. IA2I Here, the profile from the lower side 
of a sub-band exhibits a significantly suppressed first com¬ 
ponent, while that from the upper side is affected mostly at 
the phase of the trailing component. 
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Figure A2. Pulse profiles formed with a bandwidth of 3 MHz 
at the lower (left) and upper (right) side of the 1285 MHz sub¬ 
band, respectively, from the data shown in Fig. [ATI The central 
frequencies are 1277 and 1293 MHz, respectively. 
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Figure A3. Polarisation profile of the data in Fig. lAlj averaging 
over all frequencies. 


Fig. lA3l shows the polarisation profile of data in Fig. IaTI 
It can be seen that the linear and circular components in 
general retain their original shape, while the total intensity 
profile is greatly suppressed and lower than the linear com¬ 
ponent at most of the on-pulse phases including the valley 
region between the two main peaks. This explains why the 
WSRT profile shown in Fig. [T] is slightly over-polarised. 
























